pro test_camb_spec

  ; WMAP 5 spectral index
  ns = 0.963

  ; read Aquarius power spectrum
  readcol,'billennium_powspec.txt',Aq_k,Aq_k3Pk

  ; test power spectrum method
  readcol,'test_matterpower_wmap5.dat',CDV_k,CDV_untilted_Pk

  CDV_Pk   = CDV_k^(1.-ns)*CDV_untilted_Pk
  CDV_k3Pk = CDV_k^3*CDV_Pk 
  
  set_plot,'x'
  window,0,xs=512,ys=512
  plot,10.^Aq_k,10.^Aq_k3Pk,/xl,/yl
  oplot,CDV_k,CDV_k3Pk,linestyle=1,thick=5.0

  ; we want to extrapolate out to very long wavelength, low-k
  ; to enable converged integration when setting sigma_8, so sample
  ; gradient from log_10(k)= [-4,-3] and apply to [-5,-4]
  ; sample evenly in log_10(k) with 100 steps 
  i0 = value_locate(CDV_k,1.0e-4)
  i1 = value_locate(CDV_k,1.0e-3)
  gradient = (alog10(CDV_k3Pk[i1])-alog10(CDV_k3Pk[i0]))/(alog10(CDV_k[i1])-alog10(CDV_k[i0]))
  CDV_k_extend    = -5 + findgen(100)/100.
  delta_x         = CDV_k_extend - alog10(CDV_k[0])
  delta_y         = delta_x*gradient
  CDV_k3Pk_extend = alog10(CDV_k3Pk[0]) + delta_y  
  CDV_k_extend    = 10.^CDV_k_extend
  CDV_k3Pk_extend = 10.^CDV_k3Pk_extend

  ; turn into single dataset
  CDV_k    = alog10([CDV_k_extend,CDV_k])
  CDV_k3Pk = alog10([CDV_k3Pk_extend,CDV_k3Pk]) - 1.3

  set_plot,'x'
  window,2,xs=512,ys=512
  plot,CDV_k,CDV_k3Pk
  oplot,Aq_k,Aq_k3Pk,linestyle=1

  ; write this out, we'll use it
  openw,1,'wmap5_wide_powspec.dat'
  for i=0,n_elements(CDV_k)-1 do begin
     a = CDV_k[i]
     b = CDV_k3Pk[i]
     printf,1,a,b
  endfor
  close,1

; ; test power transfer function method
;  readcol,'test_transfer_wmap5.dat',CDVt_hik,CDVt_Delta_cdm,$
;          CDVt_Delta_b,CDVt_Delta_g,CDVt_Delta_r,CDVt_Delta_nu,CDVt_Delta_tot
;  CDVt_Tmatter = CDVt_Delta_tot
;  CDVt_untilted_Pk = CDVt_hik^ns * CDVt_Tmatter^2
;  CDVt_Pk   = CDVt_hik^(1-ns)*CDVt_untilted_Pk
;  CDVt_k3Pk = CDVt_hik^3*CDVt_Pk
;
;  set_plot,'x'
;  window,1,xs=512,ys=512
;  plot,10.^Aq_k,10.^Aq_k3Pk,/xl,/yl
;  oplot,CDVt_hik,CDVt_k3Pk/2.0e9,linestyle=1
;
  stop

end
